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We present a time-dependent density-functional method able to describe the photoelectron spectrum of atoms 
and molecules when excited by laser pulses. This computationally feasible scheme is based on a geometrical 
partitioning that efficiently gives access to photoelectron spectroscopy in time-dependent density-functional 
calculations. By using a geometrical approach, we provide a simple description of momentum-resolved photoe- 
mission including multi-photon effects. The approach is validated by comparison with results in the literature 
and exact calculations. Furthermore, we present numerical photoelectron angular distributions for randomly 
oriented nitrogen molecules in a short near infrared intense laser pulse and helium- (I) angular spectra for 
aligned carbon monoxide and benzene. 
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I. INTRODUCTION 

Photoelectron spectroscopy is a widely used tech- 
nique to analyze the electronic structure of complex sys- 
tems!^ The advent of intense ultra-short laser sources 
has extended the range of applicability of this technique 
to a vast variety of non-linear phenomena like high- 
harmonic generation, above-threshold ionization (ATI), 
bond softening and vibrational population trapping. 3 
Furthermore, it turned attosecond time-resolved pump- 
probe photoelectron spectroscopy into a powerful tech- 
nique for the characterization of excited-states dynam- 
ics in nano-structures and biological systems.^ Angular- 
resolved ultraviolet photoelectron spectroscopy is by 
now established as a powerful technique for study- 
ing geometrical and electronic properties of organic 
thin films Time-resolved information from streaking 
spectrograms,^ shearing interferograms, 8 photoelectron 
diffraction,^ photoelectron holography, 10 etc. hold the 
promise of wavefunction reconstruction together with 
the ability to follow the ultrafast dynamics of electronic 
wave-packets. Clearly, to complement all these experi- 
mental advances, and to help to interpret and understand 
the wealth of new data, there is the need for ab-initio the- 
ories able to provide (time-resolved) photoelectron spec- 
tra (PES) and photoelectron angular distributions (PAD) 
for increasingly complex atomic and molecular systems 
subject to arbitrary perturbations (laser intensity and 
shape). 
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Photoelectron spectroscopy is a general term which 
refers to all experimental techniques based on the pho- 
toelectric effect. In photoemission experiments a light 
beam is focused on a sample, transferring energy to the 
electrons. For low light intensities an electron can ab- 
sorb a single photon and escape from the sample with a 
maximum kinetic energy fiw — Ip (where uj is the photon 
angular frequency and Ip the first ionization potential 
of the system) while for high intensities electron dynam- 
ics can be interpreted considering a three-step modelP^l 
This model provides a semiclassical picture in terms of 
ionization followed by free electron propagation in the 
laser field with return to the parent ion, and rescattering. 
Such rescattering processes are the source of many inter- 
esting physical phenomena. In the case of long pulses, 
for instance, multiple photons can be absorbed resulting 
in emerging kinetic energies of sfvuj — Ip — Up (where s 
is the number of photons absorbed, Up = e 2 /Auo 2 is the 
ponderomotive energy and e the electric field amplitude) 
forming the so called ATI peaks in the resulting pho- 
toelectron spectrum. In all cases the observable is the 
escaping electron momentum measured at the detector. 

In general, the interaction between electrons in an 
atom or molecule and a laser field is difficult to treat 
theoretically, and several approximations are usually per- 
formed. Clearly, a full many-body description of PES is 
prohibitive, excegt for the case of few (one or two) elec- 
tron systems p2HH] As a consequence, the direct solution 
of the time dependent Schrodinger equation (TDSE) in 
the so-called single- active electron (SAE) approximation 
is a standard investigation tool for many strong-field ef- 
fects in atoms and dimers and represents the benchmark 
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for analytic and semi-analytic models P * 10 * 15 ^ ^ Pertur- 
bative approaches based on the standard Fermi golden 
rule are usually employed. For weak lasers, plane wave 
methods 5 and the independent atomic center approxima- 
tion^ nave been applied, while in the strong field regime, 
Floquet theory, the strong -field approximation^ 2 ^ and 
semiclassical method o 11 ! 29 ! 30 ! are routinely used. 

From a numerical point of view, it would be highly de- 
sirable to have a PES theory based on time-dependent 
density functional theory (TDDFT ) 31 1 321 where the com- 
plex many-body problem is described in terms of a fic- 
titious single-electron system. For a given initial many 
body state, TDDFT maps the whole many-body problem 
into the time dependence of the density from which all 
physical properties can be obtained. The method is in 
principle exact, but in practice approximations have to 
be made for the unknown exchange-correlation functional 
as well as for specific density-functionals providing phys- 
ical observables. This latter issue is much less studied 
than the former, and to the best of our knowledge a for- 
mal derivation of momentum-resolved PES from the time 
dependent density has not been performed up to now. 
In any case, several works were published addressing 
the problem of single and multiple ionization processes 
within TDDFT. For example, ionization rates were cal- 
culated for atoms and molecules j 33 " 37 ! and TDDFT with 
the sampling point method (SPM) has been em ployed in 
the study of PES and PAD for sodium clusters.^™ 

In this work, besides presenting a formal derivation of a 
photoelectron orbital functional, we report on a new and 
physically sound scheme to compute PES of interacting 
electronic systems in terms of the time-dependent single 
electron Kohn-Sham (KS) wavefunctions. The scheme re- 
lies on geometrical considerations and is based on a split- 
ting technique P^HUS The idea is based on the of partition- 
ing of space in two regions (see Fig. [I] below): In the inner 
region, the KS wave function is obtained by solving the 
TDDFT equations numerically; in the outer region, elec- 
trons are considered as free particles, the Coulomb inter- 
action is neglected, and the wavefunction is propagated 
analytically with only the laser field. Electrons flowing 
from the inner region to the outer region are recorded and 
coherently summed up to give the final result. In addition 
to the adaptation of the traditional splitting procedure 
to TDDFT, we propose a novel scheme where electrons 
can seamlessly drift from one region to the other and 
spurious reflections are greatly suppressed. This proce- 
dure allows us to reduce considerably the spatial extent 
of the simulation box without damaging the accuracy of 
the method. 

The rest of this Article is organized as follows. The 
formalism for describing photoelectrons in TDDFT is de- 
lineated in Sect. [Til In order to make contact with the 
literature, we first give a brief introduction to the state- 
of-the-art for the ab-initio calculation of PES for atomic 



case of effective single-particle theories like TDDFT in 
Sect. |IIB| In Sect. |II C| we introduce the mask method, 
an efficient propagation scheme based on space partition- 
ing. 

Three applications of the mask method are presented 
in Sect. 
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One application deals with the hydrogen 
atom and illustrates the different mask methods in a 
simple one-dimensional model also in comparison with 
the sampling point method.^ The above threshold ion- 
ization of three-dimensional hydrogen is examined and 
compared with values from the literature. In the second 
application we illustrate PADs from randomly oriented 
nitrogen molecules in a strong near-infrared ultra-short 
laser pulse. Comparison with the experiment and molec- 
ular strong- field approximation is discussed.^ The third 
application of the method regards helium- (I) (wavelength 
58 nm) PADs for oriented carbon monoxide and benzene. 
Results are discussed in comparison with the plane wave 



and molecular systems. In Sect. |II A| we introduce the 
geometrical approach in the context of quantum phase- 
space. The phase-space approach is then derived in the 



approximation. Finally, in Sect.|IV|we discuss the results 
and present the conclusions. 

All our numerical calculations were performed with the 
real-time, real-space TDDFT code Oct opus, f ree ly 
available under the GNU public license. Atomic units 
are used throughout unless otherwise indicated. 



II. MODELING PHOTOELECTRON SPECTRA 

In order to put in perspective the results of the present 
Article, we will give a brief introduction on the status 
of the principal techniques available for ab-initio PES 
calculations. We start our description with the methods 
employed to study one-electron systems. 

For one-electron systems PES can be calculated ex- 
actly from the direct solution of the TDSE. Several meth- 
ods have been employed to extract PES information from 
the solution of the TDSE. The most direct and intuitive 
way is via direct projection methods where the PES is 
obtained by projecting the wave function at the end of 
the pulse onto the eigenstates describing the continuum. 
These eigenstates are extracted through the direct di- 
agonalization of the Hamiltonian without including the 
interaction with the field. The momentum probability 
distribution can then be easily obtained from the Fourier 
transform of the continuum part of the time-dependent 
wavefunction! 2 ^ 

Another approach, that avoids the calculation of the 
full continuum spectrum, involves the analysis of the ex- 
act wavefunction |^) after the laser pulse via a resolvent 
technique! 15 1 26 1 In this case, the energy resolved PES is 
given by the direct projection on out-going wavefunc- 
tions with P{E) = \(<$>(E)\^}\ 2 = (*|J&(J57)|*>, where 
<&(E) denotes an out-going (unbound) electron of energy 
E of the laser-free Hamiltonian, and D(E) is the corre- 
sponding projection operator that can be conveniently 
approximated .E^H 

Normally, one needs accurate wave functions in a large 
space domain to obtain the correct distribution of the 
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ejected electrons. This is because the unbound parts of 
the wave packet spread out of the core region, and con- 
ventional expressions for the transition amplitude need 
these parts of the wave function. Solving the TDSE 
within all the required volume in space can easily become 
a very difficult computational problem. Several tech- 
niques were developed during the years to solve the prob- 
lem. For simple cases these difficulties can be overcome 
by the use of spherical coordinates. Geometrical splitting 
techniques have also been employedJ^^ Furthermore, 
formulations in the Kramers-Henneberger frame of refer- 
ence^ and in momentum-space^ led to calculations with 
remarkable high precision. Recently a promising surface 
flux method has also been proposed.^ 

The exact solution of the TDSE in three dimensions 
for more than two electrons is unfeasible and the limit 
rises to four electrons for one-dimensional models)^ Due 
to this limitation basically all ab-initio calculations for 
multi-electron systems are preformed under the SAE ap- 
proximation. In the SAE only one electron interacts with 
the external field while the other electrons are frozenpH 
and the TDSE is thus solved only for the active electron. 
This approximation has been successfully employed in 
several photoemis sion st udies for atoms and molecules in 
strong laser fields EEOH] However, the failure of this sim- 
ple model to describe multi-electron (correlation) effects 
calls for better schemes.^ 

The inclusion of exchange-correlation effects for a sys- 
tem of many interacting electrons can be achieved within 
TDDFT while keeping the simplicity of working with a 
set of time-dependent (fictitious) single-particle orbitals. 
In spite of transferring all the many-body problem into 
an unknown exchange-correlation functional, the lack of 
a density functional providing the electron emission prob- 
ability is a major limitation for a direct access to photo- 
electron observables from the time evolution of the den- 
sity (note that, in spite of the Runge- Gross theorem^ 
stating that all observables are functionals of the time 
dependent density, in practice we know few observables 
that can be written in terms of the time-dependent den- 
sity, one example being the absorption spectra). 

There has been some attempts to describe PES and 
multiple ionization processes with TDDFT in the stan- 
dard adiabatic approximation.^^ All these works use 
boundary absorbers to separate the bound and contin- 
uum part of the many-body wavefunction. The emission 
probability is then correlated with the time dependence 
of the number of bound electrons. 

An alternative and simple scheme is provided by the 
SPM. 38 Here the idea is to record single-particle wave- 
functions in time at a fixed sampling point rs away from 
the core. The time Fourier transform of the wavefunc- 
tion recorded at rs represents the probability of having 
an electron in rs with energy E. The probability to de- 
tect one electron with energy E in rs is then given by 



the sum over all occupied orbitals: 



P rs (E) = J2\Mrs,E)f 



(1) 



This method is easy to implement, can be extended to 
give also angular information,^ and is also clearly appli- 
cable to the TDSE in the SAE. However, it lacks formal 
derivation as it is directly based on Kohn-Sham wave- 
functions without a direct connection to the many-body 
state. Furthermore, it is strongly dependent on the po- 
sition of the sampling point and the minimum distance. 
This distance sometimes turns out to be quite large in or- 
der to avoid artifacts, and is strongly dependent on the 
laser pulse properties. We discuss further details con- 



cerning this method in Sect. Ill A 



In the following we present an alternative method in- 
spired by geometrical splitting and derive it from a phase- 
space point of view. The method can be naturally con- 
verged by increasing the size of the different simulation 
boxes. 



A. Phase-space geometrical interpretation 

An intuitive description of photoelectron experiments 
can be obtained resorting to a phase-space picture. Ex- 
perimental detectors are able to measure photoelectron 
velocity with a certain angular distribution for a sequence 
of ionization processes with similar initial conditions. 
The quantity available at the detector is therefore con- 
nected to the probability to register an electron with a 
given momentum p at a certain position r. From this 
consideration it would be tempting to interpret photoe- 
mission experiments with a joint probability distribution 
in the phase-space (r, p). Such a classical picture how- 
ever conflicts with the fundamental quantum mechan- 
ics notion of the impossibility to simultaneously measure 
momentum and position, and prevents us from proceed- 
ing in this direction. A link between the classical and 
quantum picture is needed beforehand. 

In order to make a connection to a microscopic descrip- 
tion it turns out to be convenient to extend the classi- 
cal concept of phase-space distributions to the quantum 
realm. A common prescription comes from the Wigner 
transform of the one-body density matrix with respect to 
the center of mass R = (r + r')/2 and relative s = r — r' 
coordinates. The d-dimensional (here and after d < 3) 
transform is defined as 



iu(R,p,t) 



with 



/ 



ds 

(2^1 



e ip ' s p(R + s/2,R-s/2,i), (2) 



P(r,r',£) = J dr 2 ...dr J v*(r,r2,...,r JV ,t) 

xV*(r',v 2 ,...,r N ,t), (3) 
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FIG. 1. (Color online) Schematic description of (a) the parti- 
tioning of space for the phase space method and (b) the mask 
method. Region A is the interaction region, B is the Volkov 
propagation region and C is the overlap region where and 
mix under the mask function. 



being the one-body density matrix, and 
^(ri, r2, . . . , rjvj t) the TV-body wavefunction of the 
system at time t. The Wigner function defined above 
is normalized and its integral over the whole space 
(momentum) gives the probability to find an electron 
with momentum p (position R). As the uncertainty 
principle prevents the simultaneous knowledge of po- 
sition and momentum, w(R, p) cannot be a proper 
joint distribution. Moreover it can assume negative 
values due to nonclassical dynamics. Nevertheless the 
Wigner function w(R, p) constitutes a concept close to a 
probability distribution in phase space (R, p) compatible 
with quantum mechanics. 

The quantum phase-space naturally leads to a geomet- 
rical interpretation of photoemission. One could think to 
divide the space in two regions A and B as in Fig. [I] (a), 
where region B represents the region where detectors are 
positioned and A is defined as the complement of B. In 
this picture, PES can be seen as the probability to have 
an electron with given momentum in B. It is then natu- 
ral to define the momentum-resolved photoelectron spec- 
trum as 



V(p) = lim / dRw(R,p,£) 



(4) 



where the spatial integration is carried out in region B, 
and the limit t —> oo assures that region B contains all 
photoelectron contributions. From the knowledge of the 
momentum-resolved PES [cf. Eq. Q] one can access sev- 
eral different quantities by simple integration. For in- 
stance, in three dimensions (d = 3) the energy-resolved 
PES is obtained integrating over the solid angle Q: 



P(E 



= p2/ 2 )= f 

JO 



(5) 



and the photoelectron angular distribution in the system 
reference frame is given by, 



P{E=p 2 /2,9)= P dcPV(p). 
Jo 



(6) 



In spite of giving an intuitive picture of PES, Eq. Q 
is not suited for direct numerical evaluation since it re- 
quires the knowledge of the full one-body density matrix 
in the whole space. In the next section we will make a 
contact with effective single particle theories like TDDFT 
to overcome the limitations due to the knowledge of the 
many-body wavefunction. In order to avoid integration 
over the whole space an efficient evolution scheme is pre- 
sented in Sect. Ill CI 



B. Phase space interpretation within TDDFT 

TDDFT is an effective single particle theory where 
the many-body wavefunction is described by an auxil- 
iary single Slater determina nt ^k s( t i, . . . , tn) built out 
of Kohn-Sham orbitals ^(r)P^^In order to simplify the 
notation, we drop the explicit time dependence from the 
wavefunctions and assume that the following equations 
are written in the limit t —> oo as prescribed by Eq. Q. 

Being represented by a single determinant, the one- 
body Kohn-Sham density matrix is given by 



pKs(r,r') 



occ. 



^(r)Vi(r') 



(7) 



where the sum in carried out over all occupied orbitals. 
Performing a decomposition of each orbital according to 
the partition of Fig. [I] (a) we obtain 



(8) 



where ^^(r) is the part of the wavefunction describing 
states localized in A and ^^(r) is the ionized contri- 
bution measured at the detector in B. The one-body 
density matrix can now be accordingly decomposed as a 
sum of four terms 

occ. 

PKs(r,r>) = [tPaA*WaM) +^ ) i(r)Vs ) i(r') 

i=l 

+ ^b^vWaM) + ^bA^WbM)] ■ ( 9 ) 

From Eq. (|9) we can build the KS Wigner function de- 
fined in Eq. (2k and obtain the momentum-resolved prob- 
ability distribution by inserting it into Eq. Q. We 
note that this step involves a non-trivial approximation, 
namely that the KS one-body density matrix is a good 
approximation to the fully interacting one in region B. 
This is, however, much milder than the assumption that 
the Kohn-Sham determinant is a good approximation to 
the many-body wavefunction in region B, as it is done, 
e.g., in the SPM. 

The final result is a sum of four overlap double inte- 
grals that can be simplified further. For a detailed calcu- 
lation we refer to Appendix [Aj The first overlap integral, 
containing a product of two functions localized in A [cf. 
Eq. ([9])], is zero due to the spatial integration in B. The 
two next overlap integrals, containing mixed products of 
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wavefunctions localized in A and B, can be reduced by 
increasing the size of region A. Assuming A to be large 
enough to render these terms negligible the only integral 
we are left with is the one containing functions in 
leading to 



p(p) 



L dR I 



ds 

(2tt)5 

OCC. 



i(R+o)^( R -o)- ( 10 ) 



The approximation sign ^ is a reminder for the er- 
ror committed in discarding the mixed overlap integrals. 
Since the probability of finding an ionized electron in re- 
gion A is zero for t — » oo, we can extend the integration 
over B in Eq. (10) to the whole space. Using the integral 
properties of the Wigner transform we finally obtain 
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FIG. 2. (Color online) The mask function implements a 
smooth transition from region A to region B. The functional 
shape used for actual calculations is denned in Eq. (18). 



p(p) 



\^ B ' 



(p) 



(ii) 



where V^B,i(p) is the Fourier transform of ipB,i( r ) and 
the expression is written in the limit for t — >• oo. Equa- 
tion (11) gives an intuitive formulation of momentum- 



resolved PES as a sum of the Fourier component of each 
orbital in the detector region. It is worth to note that 
Eq. (11) is not restricted to TDDFT and can be ap- 



plied to other effective single-particle formulations such 
as time-dependent Hartree-Fock and the TDSE in the 
SAE approximation. 

The numerical evaluation of the ionization probability 
from Eq. (11) requires the knowledge of the wavefunc- 
tion after the external field has been switched off. For 
ionization processes this means that one has to deal with 
simulation boxes that extend over several hundred atomic 
units and this practically constrains the method only to 
one-dimensional calculations. In the next section we will 
derive a simple scheme to overcome this limitation mak- 
ing the present scheme applicable for realistic simulations 
of molecules and nanostructures. 



C. The mask method 

In the previous sections we described a practical way 
to evaluate the momentum-resolved PES following the 
spatial partitioning of Fig. [I] (a) and how this can be 
conveniently cast in the language of TDDFT. In this sec- 
tion we take a step further in developing an efficient time 
evolution scheme by exploiting the geometry of the prob- 
lem together with some physical assumptions. 

We start by introducing a split-evolution scheme: At 
each time t we implement a spatial partitioning of Eq. ([sj) 
as following 



= M(r)^(r,t) 
^(r,t) = [l-M(r)]^(r,t) 



(12) 



where M(r) is a smooth mask function defined to be 1 
deep in the interior of region A and outside, as shown 
in Fig. [2] Such a mask function, along with the parti- 
tions A and £?, introduces a buffer region C (technically 
handled as the outermost shell of A), where ^^(r, t) and 
^B,i( r ,t) overlap [see Fig. [l](b)]. 

We can set up a propagation scheme from time t to t' 
as following 



V^(M') = [1 - M(r)]U(t',t) hMr,*) +^B,i(r,t)] 

(13) 

where U(t',t) is the time propagator associated with the 
full Hamiltonian including the external fields. Equa- 
tion ( p~3| ) defines a recursive propagation scheme com- 
pletely equivalent to a time propagation in the whole 
space A U B. 

In typical experimental setups, detectors are situated 
far away from the sample and electrons overcoming the 
ionization barrier travel a long way before being detected. 
During their journey toward the detector, and far away 
from the molecular system, they practically evolve as free 
particles driven by an external field. It seems therefore 
a waste of resources to solve the full Schrodinger equa- 
tion for the traveling electrons while their behavior can 
be described analytically. In addition, an ideal detec- 
tor placed relatively close to the molecular region would 
measure the same PES. 

From these observations we conclude that we can re- 
duce region A to the size of the interaction region and 
assume electrons in B to be well described by non- 
interacting Volkov states. Volkov states are the exact so- 
lution of the Schrodinger equation for free electrons in an 
oscillating field. They are plane-waves and are therefore 
naturally described in momentum space. In the velocity 
gauge the Volkov time propagator is formally expressed 
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by 



U v (t',t) = exp 



i 



dr- 



Mr) 



(14) 



where the time-ordering operator is omitted for brevity 
and A(t) is the vector potential. This is equivalent to the 
use of a strong-field approximation in the outer region in 
the same spirit of the Lewenstein model.^ 

In summary, the method we propose consists in solv- 
ing numerically the real-space TDDFT equations in A 
and analytically propagating the wavefunctions residing 
in B in momentum space. In this setup region C acts 
as a communication layer between functions in A and B. 
Under this prescription, and by handling 5-functions in 
momentum space, Eq. dl3|) becomes 



^A,i(r, t') = rjA,i( r i O + VB,i(*, O 
^B,*(P, O = £m(p, tf) + &,i(p, O 



with 



(15) 



(16a) 



r7B,i(r,0 =M(r) 



dp e ip,] 



Uyit'^B^t) (16b) 



(2tt)! 

x [l-M(r)}U(t',t)ip AA (v,t) (16c) 
dre -i P -. 



/ 



(2tt)* 



(16d) 



At each time step the orbital xpA,i is evolved under the 
mask function and stored in r]A,i, forcing tja,i to be lo- 
calized in A. At the same time, the components of ipA,i 
escaping from A are collected in momentum space by 
%A,i- We then add to £a,{ the contribution of the wave- 
functions already present in B at time t by summing up 
UvipB,i- In order to allow electrons to come back from 
B to A we include r]B,i in A and correct the function in 
B by removing its Fourier components [second term in 
Eq. Jl6dl)]. 



One of the advantages of Eq. (15) is that all the spa- 



tial integrals present in rjB,i(Y,t / ) and £B,i(p, £') are per- 
formed on functions localized in C. Therefore, integrals 
over the whole space are evaluated at the cost of an in- 
tegration on the much smaller buffer region C that can 
be easily evaluated by fast Fourier transform algorithms. 
Similar considerations hold for integrals in momentum 
space under the assumption that ^-functions ipB,i(p,t) 
are localized in momentum. When region A is discretized 
on a grid, in order to avoid wavefunction wrapping at the 
boundaries and preserve numerical stability, additional 
care must be taken. In our implementation, numerical 
stability is addressed by the use of non-uniform Fourier 
transforms (see details in Appendix Bj. 



There are situations were the electron flow from B to 
A is negligible. This is the case, for instance, when A is 
large enough to contain the whole wavefunctions at the 
time when the external field has been switched off. A 
propagation at later times will see photoelectrons flow- 
ing mainly from A to B. In this situation, rjB,i and the 
corresponding correction term in £b,i can be discarded. 
The evolution scheme of Eq. (15) is thus simplified and 
becomes 



^A,;(r,£') = rjA,i(r,t f ) 



(17) 



In the folowing we will refer to Eq. (17) as the "mask 
method" (MM), and to Eq. (fl5]) as the "full mask 



method" (FMM). We note here again that, being single- 
particle propagations schemes, both MM and FMM are 
not restricted to TDDFT and can be applied to other 
effective single-particle theories. As a matter of fact an 
approach similar to Eq. (17) has already been employed 



in the propagation of the TDSE equations for atomic sys- 
tems) 16 ' 17 ' 19 ' and in TDDFT for one-dimensional models 
of metal surfaces.^ We also note that the implemen- 
tation of absorb ing b oundaries trough a mask function 
as done in Eq. (16a) can be cast in terms of an addi- 
tional imaginary potential (exterior complex scaling) in 
the Schrodinger equation. Such approach is commonly 
used in quantum optics. 

Within the MM, the evolution in A is completely un- 
affected by the wavefunctions in B and ionized electrons 
are treated uniquely in momentum space. Compared 
with the FMM, the MM is numerically more stable as 
it is not affected by boundary wrapping. In order to 
achieve the conditions where Eq. ( [l7| ) is valid may re- 
quire, however, large simulation boxes. Moreover as the 
mask function never absorbs perfectly the electrons, spu- 
rious reflections may appear. Suppression of such arti- 
facts requires a further enlargement of the buffer region. 
With FMM spurious reflections are almost negligible. 

Choosing between MM and FMM implies a tradeoff 
between computational complexity and numerical stabil- 
ity that strongly depends on the ionization dynamics of 
the process under study. In what follows, we will illus- 
trate the differences and devise a prescription to help the 
choice of the most suitable method in each specific case. 



III. APPLICATIONS 

In this section we present a few numerical applications 
of the schemes previously derived. In all calculations the 
boundary between A and B regions is chosen as a d- 
dimensional sphere implemented by the mask function: 



M (r) = < 



sin 



fjr 
\2(F 



-Rc) 



(Ra-Rc 



if r < R c 

^) iiRc<r<R A , 
if r > R A 

(18) 
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as shown in Fig. [2] Note that numerical studies (not 
presented here) revealed a weak dependence of the final 
results on the functional shape of the mask. 

The time propagation of the orbitals in A is performed 
with the enforced time-reversal symmetry evolution op- 
erator^ 



U(t+At, t) = exp 



-i—H(t + At) ] exp 



(19) 

where H is the full KS Hamiltonian, and the coupling 
with the external field is expressed in the velocity gauge. 

The first system we will study is hydrogen. In spite of 
being a one-electron system, it is a seemingly trivial case 
that has been and still is under thorough theoretical in- 
vestigation.^H^EHMl clearly we do not need TDDFT to 
study hydrogen and our numerical results are obtained 
by propagating the wavefunction with a non-interacting 
Hamiltonian. The interest in this case is focused on the 
numerical performance of different mask methods as hy- 
drogen provides a useful benchmark. 

The full TDDFT calculations performed for molecular 
nitrogen, carbon monoxide, and benzene are later pre- 
sented in Sect. |IIIB] and Sect. |III C] respectively. In these 
cases, norm-conserving Troullier-Martins pseudopoten- 
tials and the exchange-correlation LB94 potential 53 (that 
has the correct asymptotic limit for molecular systems) 
are employed. Finally, in all calculations the start- 
ing electronic structure of the molecules is calculated in 
the Born-Oppenheimer approximation at the experimen- 
tal equilibrium geometry and the time evolution is per- 
formed with fixed ions. 



A. Photoelectron spectrum of hydrogen 

As first example we study multi-photon ionization 
of a one-dimensional soft-core hydrogen atom, initially 
in the ground state, and exposed to a A = 532 nm 
(uj = 0.0856 a.u.) linearly polarized laser pulse with peak 
intensity / = 1.38 x 10 13 W/cm 2 , of the form 



A(t) = A f(t) cos(ujt) 



(20) 



where f(t) is a trapezoidal envelope function of 14 optical 
cycles with two-cycle linear ramps, constant for 10 cycles, 
and with Aq = 31.7 a.u. Here A(t) is the vector potential 
in units of the speed of light c. A soft- Coulomb potential 
V(x) = —l/y/2 J rx 2 is employed to model the electron- 
ion interaction. We propagate the electronic wavefunc- 
tion in time and then compare the energy-resolved ioniza- 
tion probability obtained from different schemes. Along 
with MM, FMM, and SPM we present results for di- 
rect evaluation of PES from Eq. In this method 
the spectrum is obtained by directly Fourier transform- 
ing the wavefunction in region B. Since the analysis is 
conducted without perturbing the evolution of the wave- 
function we will refer to it as the "passive method" (PM). 




4 8 12 

t [optical cycles] 

FIG. 3. (Color online) Evolution of the electronic density as 
a function of time p(x,t) = \ip(x,t)\ 2 for the one-dimensional 
soft-core hydrogen model. The laser pulse has angular fre- 
quency uj = 0.0856 a.u., intensity / = 1.38 x 10 13 W/cm 2 , 
and a trapezoidal envelope with 2 optical cycle linear ramp 
(one optical cycle = 1.774 fs) and 10 cycles constant center. 



This method requires the knowledge of the whole wave- 
function after the pulse has been switched off, and since 
a considerable part of the wave-packet is far away from 
the core (for the present case a box of 500 a.u. radius 
is needed for 18 optical cycles), it is viable only for one- 
dimensional calculations. Nevertheless it is important as 
it constitutes the limiting case for both MM and FMM. 

In Fig. [3j a color plot of the evolution of the electronic 
density as a function of time is shown. The electronic 
wavefunction splits into sub-packets generated at each 
laser cycle (one optical cycle = 2tt/uj = 1.774 fs). These 
wavepackets evolve in bundles and their slope correspond 
to a certain average momentum. ATI peaks are then 
formed by the build up of interfering wavepackets peri- 
odically emitted in the laser field and leading to a given 
final momentumP 

From Fig. [3] is it possible to see that electrons may be 
considered as escaped "already" at 30 a.u. away from 
the center. We set therefore Ra — 30 a.u. and calcu- 
late energy-resolved PES with the PM. As we can see 
from Fig. [4] the spectrum presents several peaks at in- 
teger multiples of uj following E — suj — Ip — Up with 
Up = Aq/4c 2 = 0.0133 a.u. being the ponderomotive 
energy, Ip = 0.5 a.u. the ionization potential, and s 
the number of absorbed photons. In this case the min- 
imum number of photons needed to exceed the ioniza- 
tion threshold is s = 6. Of course, the spectrum is only 
in qualitative agreement with three-dimensional calcu- 
lations^ as expected from a one-dimensional soft-core 
modelPHSI] 

PES calculated from MM, FMM, and SPM all agree 
as reported in Fig. [4j Numerical calculations were per- 
formed until convergence was achieved, leading to a grid 
with spacing AR = 0.4 a.u. and box sizes depending on 
the method. For MM we employed a simulation box of 
Ra = 70 a.u. and set the buffer region at Rc = 30 a.u. 
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E [a.u.] 



FIG. 4. (Color online) Energy resolved photoelectron prob- 
ability P(E) (logarithmic scale) calculated with different ap- 
proaches. The spectra are shifted by multiplying a constant 
factor for easy comparison. From bottom to t op: passive 
method Eq. |5| in green, mask method Eq. (17) in red, full 
mask method Eq. (15) in blue, and sampling point method 



Eq. |T]) in purple. Laser parameters are the same as in Fig. [3] 



In order to have energy resolution comparable with PM 
we used padding factors (see Appendix [B| P = P/v = 4 
and the total simulation time was T = 18 optical cy- 
cles. For FMM a smaller box of Ra = 40 a.u. with 
Rc = 30 a.u. is needed to converge results, and P = 8, 
P/v = 2 were needed to preserve numerical stability for 
T = 18 optical cycles. For SPM two sampling points at 
f's = —500, 500 a.u. were needed to get converge d resu lts 
with a box of 550 a.u., and a complex absorber^^ at 
49 a.u. from the boundaries of the box. In addition, a to- 
tal time of T — 74 optical cycles was required to collect all 
the wave packets. The need for such a huge box resides on 
the working conditions of SPM. In order to avoid spuri- 
ous effects, the sampling points must be set at a distance 
such that the density front arrives after the external field 
has been switched off. Therefore the longer the pulse the 
further away the sampling points must be set. For these 
laser parameters one could rank each method according 
to increasing numerical cost starting from MM, followed 
by FMM, PM, and SPM. 

As a second example we study the ionization of this 
one dimensional hydrogen atom by an ultra-short intense 
infrared laser. We employ a single two-cycle pulse of 
wavelength A = 800 nm (uj = 0.057 a.u.), intensity / = 
2.5 x 10 14 W/cm 2 , and envelope 



/(*) 



- sin(c^/27V c ) 2 if < t < 2ttN c /uj 
if t > 2ttN c /uj 



(21) 



with N c = 2 and A = 225.8 a.u. 

Due to the laser strength and long wavelength, the 
electron evolution shown Fig. [5] is quite different from 




0.01 



1 2 3 

t [optical cycles] 



^0 



FIG. 5. (Color online) Density evolution p(x,t) for a one- 
dimensional soft-core hydrogen and a two-cycle sin 2 laser 
pulse with angular frequency uj = 0.057 a.u., and intensity 
I = 2.5 X 10 14 W/cm 2 . Here one optical cycle = 2.66 fs. 




FIG. 6. (Color online) Energy-resolved photoelectron P(E) 
yield from different approaches. Spectra are shifted by a con- 
stant factor. Order and color coding is the same as in Fig. [4] 
Laser parameters are described in the caption of Fig. [5] 



the one presented before. Electrons ejected from the core 
are driven by the laser and follow wide trajectories be- 
fore returning to the parent ion. Such trajectories can be 
understood in the context of the semiclassical modeJ^ 
where released electrons move as a free particle in a time- 
dependent field with a maximum oscillation amplitude 
of xq = 2Aq/ujc = 57.8 a.u. Electrons ejected near a 
maximum of the electric field e(t) = —dA(i)/dt are the 
ones gaining the most kinetic energy and are therefore 
responsible for the fast emerging electrons after rescat- 
tering with the core. 

In Fig. [6] we show the energy-resolved PES for different 
methods. Here the spectra appear to be very far from any 
ATI structure due to short duration of the laser pulse and 
is characterized by some irregular maxima and minimaP^ 
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t [optical cycles] 

FIG. 8. (Color online) Photoelectron angular distribution 
FIG. 7. (Color online) FMM density evolution |^a(x,£)| 2 of P(E,0) (logarithmic scale) of hydrogen for a 20 cycle sin 2 
a one-dimensional soft-core hydrogen. Laser parameters are laser pulse of wavelength A = 800 nm, and intensity I = 
the same as in Fig. [5] Dashed lines indicate the edges of the 5 x 10 13 W/cm 2 polarized along the x axis, 
simulation box. 



The characteristic features of the ionization dynamics is 
strongly dependent on the detailed shape of the pulse as 
one can easily imagine by inspecting the asymmetry in 
the electron ejection from Fig. [5] Due to these dynamics, 
a dramatic carrier envelope phase dependence for such 
short pulses is expected. 

All the different methods result in similar spectra but 
with different parameters. In PM we set Ra = 50 a.u. 
and a box of radius R = 700 a.u. is needed to contain 
the wave function after T = 4 optical cycles (one optical 
cycle = 2.66 fs). For MM R A = 200 a.u., R c = 40 a.u., 
and the padding factors are P = 2, P/v = 4. Here the 
value for Ra is dictated by the width of the buffer region 
which needs to be wide enough to prevent spurious re- 
flections. A considerably smaller box is needed for FMM, 
where Ra = 60 a.u., Rc = 40 a.u., P = 4, and P/v = 2. 
In this case one can reconstruct the total density in A 
by evaluating |^(x,t)| 2 via Eq. (15) and compare it to 
the exact evolution. As one can see in Fig. [7] the recon- 
structed density displays a behavior remarkably similar 
to the exact one of Fig. [5] but with a considerably reduced 
computational cost. SMP requires sampling points at 
rs = —130, 130 a.u. in a box of radius R = 200 a.u. with 
49 a.u. wide complex adsorbers, and for a total time of 
T = 7 optical cycles. 

The possibility to use relatively small simulation boxes 
is especially important for three-dimensional calculations 
where the computational cost scales with the third power 
of the box size. Both mask methods are practicable op- 
tions for 3D simulations and the advantage of using FMM 
with respect to MM is driven by the electron dynamics. 
For long laser pulses MM appears to be more stable and 
it is a better choice than FMM, while for short pulses 
with large electron oscillations FMM can be more per- 
formant. SPM is a viable option for short pulses and 
small values for the oscillations. 

As a last example, we present ATI of a real three- 
dimensional hydrogen atom subject to a long infrared 



pulse. We employ a laser linearly polarized along the 
x-axis with wavelength A = 800 nm, intensity J = 5 x 
10 13 W/cm 2 , pulse shape of the form ^ with N c = 20 
and Ao = 91.3 a.u. Due to to the pulse length, MM 
appears to be the most appropriate choice in this case. 
In the calculation Ra = 60 a.u., Rc = 50 a.u., P = 1 
and P/v = 8. 

In Fig. m we show a high-resolution density plot of the 
PAD P(E,6) defined in Eq. The radial distance de- 
notes the photoelectron energy while the angle indicates 
the direction of emission with respect to the laser polar- 
ization. The color density is plotted in logarithmic scale 
and represents the values of P(E,0). 

The photoelectron energy- angular distribution dis- 
plays complex interference patterns. The pattern shape 
compares favorab ly with similar calculations in the lit- 
erature P 4 * 25 * 44 * 59 ! It consists of a series of rings with fine 
structures. Each ring represents the angular distribution 
of the photoelectron ATI peaks. The spacing of adjacent 
rings equals the photon energy uj = 0.057 a.u. Photo- 
electrons are emitted mainly along the laser polarization, 
and the left-right symmetry of the rings indicates that 
the photoelectrons do not present any preferential ejec- 
tion side with respect to the polarization axis. The first 
ring corresponds to the angular distribution of the first 
ATI peak. It presents a peculiar nodal pattern that is 
induced by the long-range Coulomb potential and is re- 
lated to the fact that the ATI peak is determined by one 
dominant partial wave in the final stated The number of 
the stripes equals the angular momentum quantum num- 
ber of the dominant partial wave in the final state plus 
one Jlo] j n |gj first ring contains six stripes and 
the dominant final state has angular momentum quan- 
tum number of 5. The pattern of the energy- angular 
distribution and the stripe number of the first ring are in 
good agreement with those in the literature! 24 * 44 * As for 
the fine structures, we observe that while the main ring 
pattern is already formed in the first half of the pulse, 
the fine structure builds up until the end of the pulse. 
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This supports the hypothesis that such structures are in- 
duced by the coherence of the two contributions from the 
leading and trailing edges of the pulse envelope!^ 



B. N2 under a few-cycle infrared laser pulse 

In this section we compare theoretical and experimen- 
tal angular resolved photoelectron probabilities for ran- 
domly oriented N2 molecules. We choose the laser pa- 
rameters according to experiment, 28 i.e., we employ a 
N c = 6 cycle pulse of wavelength A = 750 nm (u = 0.06 
a.u.), intensity / = 4.3 x 10 13 W/cm 2 . A laser shape 



A(t) 



cos(#))sin(wt) 



if < t < 2ttN c /lu 
if t > 2irN c /u 

(22) 

for the vector potential should lead to an electric field 
similar to the one employed in the experiment with zero 
carrier envelope phase. 

In Fig. M (a) the experimental photoelectron probabil- 
ity P(E : 0) is plotted in logarithmic scale as a function 
of the energy and the angle with respect to the laser po- 
larization in the laboratory frame. Electrons are mainly 
emitted at small angles, and, due to the short nature of 
the pulse electron emission is asymmetric along the laser 
polarization axis (at angles close to 0° and 180°). 

We performed TDDFT calculations for different an- 
gles 0l between the molecular axis and the laser po- 
larization. The molecular geometry was set at the ex- 
perimental equilibrium interatomic distance R$ = 2.074 
a.u. The Kohn-Sham wavefunctions were expanded in 
real space with spacing Ar = 0.38 a.u. in a simulation 
box of Ra = 35 a.u. The photoelectron spectra were 
calculated with FMM having Rq = 25 a.u., and padding 
factors P = 1, and P/v = 4. 

In Fig. [To] the logarithmic ionization probability 
Pq l (£", 0) is plotted as a function of energy E and angle 
measured from the laser polarization axis for different 
values of 0l. As the molecular orientation decreases from 
90° < 0l < 30° we observe an increasing suppression of 
the emission together with a shift of the maximum that 
moves away from the laser polarization axis. For 0l = 0° 
the emission is highly enhanced for all angles and peaked 
along the laser direction. The signature of multi-center 
emission interference has been predicted to be particu- 
larly marked when the l aser p olarization is perpendicu- 
lar to the molecular axi J 61 * 62 1 (i.e. 0l = 90°). However, 
the lowest point in energy of such a pattern is predicted 
for = 90° and E = n 2 /2Rl w 31 eV, way above the 
energy window of observable photoelectrons produced by 
our laser. A stronger and longer laser pulse would be re- 
quired to extend the rescattering plateau toward higher 
energies and therefore to reveal the pattern.^ 

In order to reproduce the experimental P(E, 0), an av- 
erage over all the possible molecular orientations should 
be performed. Due to the axial symmetry of the molecule 
we can restrict the average to < 0l < 90° and integrate 
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FIG. 9. (Color online) Photoelectron angle and energy re- 
solved probability P(E,9) (log scale) in the laboratory frame 
for randomly oriented N2 molecules in a 6 cycles infrared 
laser pulse with A = 750 nm, with intensity / = 4.3 x 10 13 
W/crn . The angle 6 is measured from the laser polarization 
axis. The different panels represent P(E, 0) (spanning 3.4 or- 
ders of magnitude), from (a) experiment, (b) calculated with 
TDDFT and FMM, and (c) calculated with modified molecu- 
lar strong field approximation. Panels (a) and (c) are adapted 
from Ref . M 



all the contributions with the proper probability weig ht 41 



^.90° 
^0 



d9 L sm6 L Pe L (E,6). 



(23) 



We evaluate Eq. (23) by discretizing the integral in a 
sum for 6 L = 0°, 30°, 60°, 90°, and display the result in 
Fig. [9](b). Even in this crude approximation, and with- 
out taking into account focal averaging, the agreement 
with the experiment is satisfactory and compares favor- 
ably to the molecular strong field approximation shown 
in Fig.[9](c). The agreement deteriorates for low energies 
where the importance of the Coulomb tail is enhanced as 
it is not fully accounted due to limited dimensions of the 
simulation box. As a matter of fact the agreement greatly 
increases for higher energies. 
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FIG. 10. (Color online) Photoelectron angle and energy 
resolved probability Pq l (E, 0) (log scale) for aligned N 2 
molecules and different laser polarization directions 0l- (a) 
Q L = 90°, (b) L = 60°, (c) L = 30°, (d) L = 0°. Here L 
is the angle between the laser polarization direction and the 
molecular axis. Laser parameters are the same as in Fig. [5] 



C. He-(l) PADs for carbon monoxide and benzene 

In this section we deal with UV (uj = 0.78 a.u.) an- 
gular resolved photoemission triggered by weak lasers. 
When the external field is weak, non-linear effects can 
be discarded and first order perturbation theory can be 
applied. In this situation, the momentum resolved PES 
can be evaluated by Fermi's golden rule as 

P(p) oc K*/l A o • P\*i)\ 2 5(E f -Ei-u,), (24) 

i 

where |\E^) (|^/)) is the initial (final) many-body wave- 
function of the system and Ao is the laser polarization 
axis. The difficulty in evaluating Eq. (24) lies in the 
proper treatment of the final state, which in principle be- 
longs to the continuum of the same Hamiltonian of 
In the simplest approach, it is approximated by a plane 
wave (PW) . In this approximation the square root of the 
momentum-resolved PES is proportional to the sum of 
the Fourier transforms of the initial state wavefunctions 
^i(p) corrected by a geometrical factor |A • p| 



y/Pjpj OC |Ao • Pi X |#;(p)|. 



(25) 



If photoemission peaks are well resolved in momentum, 
individual initial states can be selectively measured. In 
this case a correspondence between momentum-resolved 
PES and electronic states in reciprocal space can be es- 
tablished. The range of applicability of the PW ap- 
proximation has been discussed in the literatureP It has 




FIG. 11. Photoemission geometries for oriented (a) benzene 
and (b) CO molecules. 



been postulated that Eq. p5[ ) should be valid for (i) 
7r-conjugated planar molecules, (ii) constituted by light 
atoms (H, C, N, O) and for (hi) photoelectrons emerg- 
ing with momentum p almost parallel to the polarization 
axis. 

Here we restrict ourselves to photoemission from the 
highest occupied molecular orbital (HOMO). In this case 
Eq. (25) becomes 



VP* (p) cx |A • p| x |# H (p)|, 



(26) 



the subscript H indicating HOMO-related quantities. 
We compare ab-initio TDDFT and PW PADs evaluated 
at fixed momentum \ph\ = V^Eh with Eh = uj — Eb 
being the kinetic energy of photoelectrons emitted from 
the HOMO and Eb its binding energy. 

TDDFT numerical calculations are carried out on a 
grid with spacing Ar = 0.28 a.u. for benzene and 
Ar = 0.38 a.u. for CO, in a simulation box of Ra = 
30 a.u.. Photoelectron spectra are calculated using MM 
with Rc = 20 a.u. and padding factors P = 1, P/v = 8. 
A 40 cycles pulse with 8 cycle ramp at the He- (I) fre- 
quency uj = 0.78 a.u. and intensity / = 1 x 10 8 W/cm 2 
is employed. 

We begin presenting the case of benzene since it con- 
stitutes the smallest molecule meeting all the conditions 
for Eq. (26) to be valid. Results for molecules oriented 
according to Fig. [IT] (a), evaluated at Eh = 0.363 a.u., 
and two different laser polarizations Aq = a x , a 2 with 
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ai = (1, 0, 0), a 2 = l/\/3 x (1, 1, 1), are shown in Fig. 
In the case where the laser is polarized along the x axis 
[see Fig. 12 (b)], PAD presents a four lobes symmetry sep- 



arated by three horizontal and two vertical nodal lines. 
This structure is reminiscent of the HOMO 7r-symmetry 
with the nodal line at = 90° corresponding to the nodes 
of the orbital on the x-y plane. Information on the orien- 
tation of the molecular plane could then be inferred from 
the inspection of this nodal line in the PAD. A similar 
feature can be observed also in the case of an off-plane 
polarization as shown in Fig. [l2] (d). In this case, how- 
ever, the laser can also excite a-orbitals and the nodal 
line at 6 = 90° is partially washed out. The other nodal 
lines can be understood in term of zeros of the polar- 
ization factor |Aq • p| and are thus purely geometrical. 
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FIG. 12. (Color online) He-(I) PADs for aligned benzene 
molecules. We compare PADs from PW |A • p||^ff (p)| (left 
column) and TDDFT y/ Ph (p) (right column) on a sphere at 
constant kinetic energy Eh — 0.363 a.u. for different laser po- 
larizations Ao (see text for details) . Values on the sphere are 
normalized to unity. We used a 40 cycles (8 cycles ramp) UV 
trapezoidal laser pulse with A = 58 nm (cj = 0.78 a.u.), and in- 
tensity I = 1x10 s W/cm 2 . In the top row A = ai = (1,0,0), 
and in the bottom row Ao = a.2 = 1/V3x (1,1,1). White tics 
indicate the intersection of the laser polarization axis with the 
sphere at constant kinetic energy Eh- T he g eometry of the 
photoemission process is indicated in Fig. 11 (a). 



A PW approximation of the photoelectron distribution 
given by Eq. (26) qualitatively reproduce the ab-initio 
results as shown in Fig. [l2] (a) and (d). According to 
condition (iii) a quantitative agreement is reached only 
for directions parallel to the polarization axis. 

A different behavior is expected in the case of CO. 
Photoelectrons with kinetic energy of Eh = 0.261 a.u. 
are show in Fig. 13 In this case, condition (i) (i.e. ir- 



conjugated molecule) is not fulfilled and a worse agree- 
ment between ab-initio and PW calculations is expected. 
The quality of the agreement can be assessed by com- 
paring the left and right columns of Fig. [13] Here, the 
weak angular variation of |\j/#(p)| is completely masked 
by the polarization factor |Aq • p| [cf. Fig. 13 (a) and 
(c)]. For this reason no information on the molecular 
configuration can be recovered from a PW model. 

The situation is qualitatively different for TDDFT as, 
in this case, single atom electron emitters are fully ac- 
counted for. Here the nodal pattern is mainly gov- 
erned by the polarization factor, but, however, finger- 
prints of the molecule electronic configuration can be de- 
tected. For instance, when the laser is polarized along 
the molecular axis, an asymmetry of the photoemission 
maxima can be observed for directions parallel to ai [see 




FIG. 13. 
molecules 
as in Fig 
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(Color online) He- (I) PADs for aligned CO 
Panel ordering and laser parameters are the same 

"(b) 
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The molecule is oriented according to Fig. 
and the photoelectron spectra were evaluated on a sphere at 
E h = 0.261 a.u. 



Fig. 13 (b)]. Here the global maximum is peaked around 
(</>, 0) = (180°, 90°) corresponding to the side of the car- 
bon atom on the molecular axis [cf. Fig. 11 (b)]. These 
features can be again understood in terms of the shape 
of the HOMO. For CO, in fact, the HOMO is a a orbital 
with the electronic charge unevenly accumulated around 
the carbon atom. It is therefore natural to expect pho- 
toelectrons to be ejected mainly around the molecular 
axis and with higher probability form the side of the car- 
bon atom. This asymmetry is therefore a property of 
the electronic configuration of the molecule and gives in- 
formation about the molecular orientation itself. This 
behavior appears to be stable upon molecule rotation as 
can be observed in the case where the polarization is 
tilted with respect to the molecular axis [Ao = ai, see 
Fig. 13 (d)]. Even here the nodal structure is mainly 
dictated by the polarization factor. 



IV. CONCLUSIONS 

In this work we studied the problem of photoemission 
in finite systems with TDDFT. We presented a formal 
derivation of a photoelectron density functional from a 
phase-space approach to photoemission. Such a func- 
tional can be directly applied to other theories based on 
a single Slater determinant and the derivation could serve 
as a base for extensions to more refined models. 

We proposed a mixed real- and momentum-space evo- 
lution scheme based on geometrical splitting. In its com- 
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plete form it allows particles to seamlessly pass back and 
forth from a real-space description to a momentum- space 
description. The ordinary splitting scheme turns out to 
be a special case of this more general method. Further- 
more, we illustrated applications of the method on four 
physical systems: hydrogen, molecular nitrogen, carbon 
monoxide and benzene. 

For hydrogen we presented a comparison of the differ- 
ent methods. We studied ATI peak formation in a one- 
dimensional model and ATI angular distributions for a 
three-dimensional case. The results turned out to be in 
good agreement with the literature. From the compari- 
son, we derived a prescription to choose the best method 
based on a classification of the electron dynamics induced 
by the external field. 

We investigated angular-resolved photoemission for 
randomly oriented N2 molecules in a short intense IR 
laser pulse. We illustrated the results for four different 
molecular orientations with respect to the laser polariza- 
tion. Owing to the symmetry of the problem we were able 
to combine the results to account for the random orien- 
tation. The spectrum for randomly oriented molecules is 
in good agreement with experimental measurements and 
is much better than the widely used strong field approx- 
imations (with one active electron) PHI 

We also studied UV angular resolved photoelectron 
spectra for oriented carbon monoxide and benzene 
molecules. We presented numerical calculations for two 
different directions of the laser polarization and com- 
pared with the plane-wave approximation. We found 
that the plane-wave approximation provides a good de- 
scription for benzene while failing for CO. Furthermore, 
we found evidence that the photoelectron angular distri- 
bution carries important information on molecular orien- 
tation. 

The successful implementation of photoelectron den- 
sity functional presented in this Article paves the road 
for interesting applications to many different systems 
for a wide range of laser parameters. To name a few, 
TDDFT PAD could provide a theoretical tool superior 
to the plane- wave and the independent atomic center ap- 
proximations to retrieve molecular adsorption orientation 
information from experiments. Atto-second pump probe 
experiments could be simulated ab-initio accounting for 
many-body effects but with great computational advan- 
tage with respect to full many-body methods and better 
physical description than SAE pictures. 
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Appendix A: Overlap integrals 

In this section we describe the details of the inclusion of 
the Kohn-Sham one-body density matrix ([£) into Eq. Q. 
The momentum- resolved photoelectron probability is the 
sum over all the occupied orbitals of four overlap integrals 

7 

occ. 

P(P) = ^27A,A,i(p)^lA,B,i(p)^^B,A,i(p)^^B,B,i(p) ■ 

2=1 

(Al) 

In order to simplify the notation we drop the orbital in- 
dex i in the overlap integrals and indicate with v = vir 
the vector v of modulus v and direction v. In addition, 
we will consider the simple case where the boundary sur- 
face between region A and B is a d-dimensional sphere 
of radius Ra- 

We start by considering the mixed overlap 

7AB(p) = j B dR I ^| eiPS ^( R +|)^(R-|) 

' B (A2) 



where the integration in B is for R > Ra [cf. Fig. [T](a)]. 
It is convenient to work in the coordinates v = 2R and 
r = R + s/2, where the integral takes the form 

/ dv / 7^ e^^-^MVMv - r) . (A3) 
Jv>2R A J {2n)2 

We substitute ipg with its Fourier integral representation 

_dk 

(27 

and after few simple steps we obtain 

(27TJ 2 Jv>2R A 

(A5) 

where we successfully disentangled the integration over 
v in the second integral. The integral on v > 2Ra can 
be rewritten as an integral over the whole space, which 
yields a d-dimensional Dirac delta, minus an integral on 
v < 2R A : 



%(u-r)= [ ^ c *-< u -^|,(k) (A4) 
J (27r)2 



dv e 



i(k+p)-v 



v>2R A 



( ^(k + p)-(M A )i ^^ |k : p| 



(A6) 
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where J n (k) is a Bessel function of the first kind. The 
second term in (A6) is a function centered in — p and 
strongly peaked in the region w = Cd/RA with C\ = 7r, 
C2 ~ 3.83, Cs ~ 4.49 being the first zeros of the Bessel 
function If the region w i s sm all enough we 

can consider the integrand in k of ( |A5[ ) constant and 
factor out of the integrand ^(— 2p — k)^(k) evaluated 
at k = — p. It is easy to see that 



dk 

(2^ 



(2R A ) 



, J d/2 (2i? A |k + p| 
|k + p|f 



1 



(A7) 



and, by plugging ( A6) in ( A5), we have that Ja,b(p) ~ 0. 
By the same reasoning we should expect jb, a (p) ~ 0. 

We now turn to the terms containing wavefunction on 
the same region. In (v, r) coordinates 



7a,a(p) 



Jv>2R A J 



dr 

(2^1 



ip-(2r-v) 



(A8) 

The product of functions localized in A is not negligible 
only for r < Ra and |v — r| < Ra- Since the integral is 
carried out for v > 2Ra we can bound | v — r| from below 
with R a \2v-t/R a \ > Ra- This leads to R A < |v-r| < 
Ra which is satisfied only on the boundary of A. Being 
a set of negligible measure we have ja, a (p) = 0. 
Once again, in (v, r) coordinates 



Jb,b(p) = I dv / -i^ e ^( 2r - v )^(r)^(v-r) 

Jv>2R A J (27TJ2 

(A9) 



u>2R A 

can be written as 



7b,b(p) = \^b(p)\ 2 - 
dr 



dv 



v<2R A 



(2tt)1 



jp-(2r-v) 



V> B (r)^(v-r) (AlO) 



where the first integration is in region A. Using the lo- 
calization of ipB we see that the integral is non-zero only 
for r > Ra and |v — r| > R A . As the integration is for 
v < 2Ra we have that Ra > |v — r| > Ra and therefore 



the double integral in Eq. (AlO) is zero. 



Appendix B: Numerical stability and Fourier integrals 



A real- space implementation of Eq. ( 15 ) involves the 
evaluation of several Fourier integrals. Such integrals 
are necessarily substituted by their discrete equivalent, 
and therefore discrete Fourier transforms (FT) and fast 
Fourier transforms (FFTs) are called into play. How- 
ever, evolution methods based on the discrete FT nat- 
urally impose periodic boundary conditions. While this 
is not presenting any particular issue for MM where FT 
are only used to map real-space wavefunctions to mo- 
mentum space, it is a source of numerical instability for 
FMM where the wavefunctions are reintroduced in the 
simulation box. 



The problem is well illustrated by the following one- 
dimensional example. Imagine a wavepacket freely prop- 
agating to an edge of the simulation box with a certain ve- 
locity. In MM, when passing trough the buffer region, the 
packet is converted by discrete FT in momentum space 
and then analytically evolved as a free particle through 
the edge of the box. In FMM as the wavefunction evolves 
in momentum spaces it is also transformed back to real 
space to account for possible charge returns. In this case, 
instead of just disappearing from one edge, by virtue of 
the discrete FT periodic boundary conditions, the same 
wavepacket will appear from the opposite side. It can be 
easily understood how such an undesirable event can cre- 
ate a feedback leading to an uncontrolled and unphysical 
build up of the density. 

This behavior can be controlled by the use of zero 
padding. As we know, the Fourier integrals in Eq. (15) 



involves functions that are, by construction, zero outside 
the buffer region C. We can therefore enlarge the inte- 
gration domain (having radius Ra) by a padding factor 
P, set the integrand to zero in the extended points, ob- 
taining the same result. As a consequence, a wavepacket 
propagating toward a boundary edge will have to run an 
enlarged virtual box of radius Ra = R A (2P — 1) before 
emerging from the other side. In addition, the smallest 
momentum represented Ap = 2tt/PRa = Ap/P in the 
discretized ipB,i(p, t') is reduced by a factor 1/P while the 
highest momentum Pmax = 7r/Ar remains unchanged. 
The price to pay here is an increased memory require- 
ment by a factor P d (where d is the dimension of the 
simulation box) and is too high for three-dimensional cal- 
culations. 

A possible way to find a better scaling is offered by the 
use of non-uniform discrete Fourier Transform and com- 
panion fast algorithm NFFTpHH NFFT allow for the 
possibility to perform Fourier integrals on unstructured 
sampling points with, for fixed accuracy, the same arith- 
metical complexity as FFT. For a detailed description 
of the algorithm we refer to the literature.^ The idea is 
to use the flexibility of NFFT to perform zero padding 
in a convenient way. Instead of allocating an enlarged 
box filled with zeros at equally spaced sample positions, 
we set only one point at RaPn (here P/v is the NFFT 
padding factor) and evaluate the Fourier integral with 
NFFT. In this way we gain numerical stability for FMM 
as long as the wavefunctions are contained in a virtual 
box of Ra = Ra(2Pn — 1) at the price of adding a num- 
ber of points that scales as d — 1 with the dimension of 
the box. If N d is the number of grid points in the simu- 
lation box, in order to perform zero padding with NFFT 
one needs to add only 2N d ~ 1 points. 

With this procedure however, not only the smallest 
momentum Ap is reduced by a factor 1 / P/v > but also the 
highest momentum p m ax — {N/2 + l)Ap is decreased 
by the same amount. This turns out to be the limit- 
ing factor in the use of NFFT to preserve numerical sta- 
bility with FMM as the enlargement factor P/v has an 
upper bound that depends on the escaping electron dy- 
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nami cs. I n fact, when we evaluate the back-action term 
Eq. (16b), we assume ^^(p, £') to be localized in mo- 



mentum and, in order to preserve numerical consistency, 
P/v must be limited by the highest momentum contained 
in i/jB,i- A combination of ordinary padding and NFFT 
padding helps to balance the tradeoff between memory 
occupancy and numerical stability. 

Finally, in MM zero padding can be used to increase 
resolution in momentum. 
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